Density Matrix Approach to Local Hilbert Space Reduction 
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We present a density matrix approach for treating systems with a large or infinite number of degrees 
of freedom per site with exact diagonalization or the density matrix renormalization group. The 
method is demonstrated on the ID Holstein model of electrons coupled to Einstein phonons. In this 
system, two or three optimized phonon modes per site give results as accurate as with 10-100 bare 
phonon levels per site. 



PACS Numbers: 74.2Q.Mn, 71.10.Fd, 71.10.Pm 

During the past decade there have been great strides in 
the development of numerical techniques for simulating 
strongly correlated systems. A significant limitation to 
many of these methods — for example, exact diagonaliza- 
tion using the Lanczos algorithm, or the density matrix 
renormalization group |]l| (DMRG) — is that they require 
a finite basis. In an electron-phonon lattice model, for 
example, the number of phonons is not conserved and 
the Hilbert space is infinite, even for a finite number of 
sites. Of course, the number of phonons can be artificially 
constrained, but for strongly coupled systems the num- 
ber of phonons needed for an accurate treatment may be 
quite large. This often severely constrains the size of the 
systems which may be studied. 

Here we present a technique for generating a controlled 
truncation of the Hilbert space, which allows the use of 
a very small local basis without significant loss of accu- 
racy. The local basis which is generated can be used in 
exact diagonalization, DMRG, or other approaches to al- 
low treatment of larger systems. The procedure is closely 
related to DMRG, in that the local basis is generated us- 
ing a density matrix, but it is simpler to implement. We 
illustrate the method on the Holstein model, a model of 
noninteracting electrons on a lattice coupled to phonons, 
with one Einstein oscillator on each site. For this model 
we show that with two or three optimized phonon modes 
per site, we obtain the same accuracy as with dozens of 
unoptimized phonon levels per site. For simplicity, the 
Hilbert space reduction technique is used here in conjunc- 
tion with exact diagonalization, but coupling it to other 
approaches, such as DMRG, would allow treatment of 
larger systems. 

The key idea of this approach is identical to the key 
idea of DMRG Q : in order to eliminate states from a part 
of a system without loss of accuracy, one should trans- 
form to the basis of eigenvectors of the reduced density 
matrix, and discard states with low probability. The key 
difference is that here the subsystem is a single site, or 
a handful of sites, rather than varying fractions of the 
entire system. To be specific, consider a many-body sys- 
tem divided into unit cells which we will call "sites" , such 
that the complete Hilbert space of the system is the outer 



product of the states of the sites. The number of states 
per site may be finite or infinite. We will consider here 
only systems with translational invariance, so that all 
sites are equivalent. Let i label the states of a particular 
single site, say site 1. Let j label the combined states of 
all of the rest of the sites. Then a wavefunction of the 
system can be written as 



(1) 



The density matrix for this site when the system is in the 
state \ip) is 



(2) 



Let Wa be the eigenvalues of p, and let Wq be the eigenvec- 
tors. The Wa are the probabilities of the states v^- If Wa 
is negligible, then the corresponding eigenvector Va can 
be discarded from the basis for the site, without affecting 
the state ijj. If one wishes to keep a limited number of 
states TO per site, then the best states to keep are the 
eigenstates of p with largest eigenvalues In the case 
of the Holstein model, we will show that all but a handful 
of these eigenstates have negligible probability. 

Usually the target state ip which one wants to repre- 
sent is the ground state. If one wants a site basis which 
represents several states, one can add each state into the 
density matrix 



(3) 



where the Oq are weights assigned to each target state, 
representing the importance of that state. Again, the 
optimal states to keep are the eigenstates of p. 

Unfortunately, in order to obtain the optimal states, 
we need the target state, which we usually we do not 
know-usually we want the optimal states to help get the 
target state. This problem can be circumvented in sev- 
eral ways. We illustrate these approaches in the Holstein 
model, which has as Hamiltonian 
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where and cea are electron creation and annihilation 
operators, 6j and are phonon creation and annihila- 
tion operators, % = cJ|Qi- + cJj^qj^ and i is the hopping 
integral, 7 is the electron-phonon coupling constant, and 
each oscillator has frequency uj. 

The first approach is illustrated in Fig. 1(a). Here one 
site of the system (the "big site") is allowed to have a 
large number of phonon states M, with M ~ 10 — 100. 
The rest of the sites have a much smaller number of 
phonon levels, m ~ 2 — 3. A set of Davidson or Lanczos 
exact diagonalizations are performed. In the first diago- 
nalization, all of the phonon states are "bare" : they are 
eigenstates of the single site phonon Hamiltonian, char- 
acterized by the frequency w. After each diagonalization, 
the density matrix for the phonon modes of the big site 
is diagonalized. The most probable m eigenstates are 
the new optimal phonon modes. These optimal phonon 
modes are used on all the other sites for the next diag- 
onalization. In each diagonalization, the big site always 
has a large number of phonon modes, so that it can gener- 
ate improved optimal modes for the next iteration. The 
diagonalizations are repeated until the optimal modes 
have converged. In the first diagonalization, the optimal 
modes are not very accurate because the bare phonon lev- 
els used on the other sites severely truncate the Hilbert 
space. Convergence takes only a few diagonalizations, 
however. 
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FIG. 1. Algorithms for constructing optimal bases, (a) 
Tiic big site iias the complete set of bare phonon levels, (b) 
The big site has the optimal levels plus a few bare levels. 



While this approach is very simple, the large number of 
states on the big site limits the size of the system which 
can be diagonalized. A more sophisticated approach is 
illustrated in Fig. 1(b). Here, the big site has both the 
optimal modes and two or three extra levels, rather than 
M bare levels. These extra levels allow improvement of 
the optimal basis. They are taken from the set of M 
bare levels but are explicitly orthogonalized to the cur- 
rent optimal modes. After a diagonalization including 
these levels, a new density matrix is formed and optimal 
modes found. These optimal modes can mix in a little of 
the extra levels to improve the basis. The next diagonal- 
ization uses different extra levels. One sweep consists of 
enough diagonalizations to include all M bare levels as 
extra levels. A couple of sweeps are needed to reach full 
convergence of the optimal levels. Each diagonalization 
uses as the starting wavefunction the converged wave- 
function from the last step. Therefore only two or three 
Davidson steps are needed for convergence, rather than 
dozens. 

A further improvement comes from forming the den- 
sity matrix for an entire site, including electron degrees of 
freedom, rather than just the phonon levels. This forms 
different optimal phonon modes for each of the four elec- 
tron states of the site, which reduces the number of states 
needed for a given accuracy. 
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FIG. 2. Ground state energy of a 4-site half-filled Holstein 
system with uj = t and 7 = 1.5t as a function of the number 
of phonon states kept on each site of the lattice. 

In Fig. 2 we show the ground state energy of a four 
site system as a function of the number of phonon modes, 
both for bare phonon levels and for optimal modes, where 
the TO optimal modes are allowed to be different for each 
electron state of a site. The energy in these restricted 
bases form upper bounds to the true ground state energy 
of the four site system. The improvement coming from 
using optimal modes is remarkable. Using only two op- 
timal modes, the energy is accurate to less than 0.1%, 
whereas keeping eleven bare modes (the largest number 
we could treat), the error is greater than 5%. 
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The very rapid convergence of the energy as a function 
of the number of optimal phonon modes is due to the 
very tight distribution of their probabihties w„- In two 
special cases (7 = or f = 0) it is possible to show that at 
most one optimal phonon level for each electron state of 
a site has a non-zero probability. For arbitrary couplings 
several optimal phonon levels can have a finite probabil- 
ity but only a few of them are significant. In the ground 
state of a 6-site Holstein model with w = t,we have found 
that the third and fourth highest eigenvalues (for a 
given electron state of a site) are smaller than 10"'^ and 
10~^, respectively, for any electron-phonon coupling 7. 
Therefore, we have never needed more than m = 2 — 3 
optimal phonon levels to get an accurate ground state. 

In Fig. 3(a) we show the four most probable phonon 
states expressed in the bare phonon basis for the groimd 
state of a 6-site Holstein lattice. The probability Wa and 
the occupation of the electronic site are also given for 
each state. It is clear that one needs at least ^ 20 bare 
phonon levels for an accurate treatment of this system. 

An additional improvement is often possible: the opti- 
mal levels may be transferable from a smaller system to a 
larger system. We find that in the Holstein model at half 
filling, the levels obtained from applying this procedure 
to a two site or four site system work very well for larger 
systems. Thus for the larger system, one needs to do only 
one diagonalization, and each site has only m levels. This 
would be the simplest way to incorporate this approach 
into DMRG: use a small system diagonalization approach 
to find optimal levels, and then use DMRG to treat very 
long chains. Away from half-filling, because only certain 
fillings are possible on small systems, it would probably 
be necessary to simultaneously target two or more states 
with different fillings in order for the optimal basis to be 
transferable to a large system. 

The form of the optimal phonon levels can tell us im- 
portant information about the system. In Fig. 3(b) we 
show optimal phonon wavefunctions as a function of the 
oscillator position q = + b for different electron-phonon 
couplings. Only the most probable optimized phonon 
level is shown for each possible electron state of a site. 
If the optimal states are allowed to be different for each 
electron state, we find that every optimal state is either 
an eigcnstatc of an oscillator with a shifted equilibrium 
position or a linear combination of two such eigenstates, 
to surprisingly high accuracy, with overlaps greater than 
96%. Unfortunately, we can not use this property to cal- 
culate a priori a basis of optimal states. 

The form of the most important phonon states shown 
in Fig. 3(b) can be understood qualitatively in the weak 
and strong coupling regimes. The behavior at intermedi- 
ate couplings interpolates smoothly between strong and 
weak coupling. In the weak coupling regime (7 = OM) 
optimal states are simply eigenstates of an oscillator with 
an equilibrium position g « 27/a; as predicted by a mean- 




FIG. 3. (a) The first four optimal phonon states expressed 
in the bare phonon basis for a 6-site half-filled Holstein system 
with LJ = t and 7 = 1.5t. (b) Optimal phonon wavefunctions 
as a function of the oscillator position q for different elec- 
tron-phonon couplings 7 for the same lattice. 

field approximation. Differences between the optimal 
phonon states for the different electronic states of a site 
arc small compared to their widths, which are determined 
by phonon fluctuations. Therefore, electron and phonons 
are almost independent and the lattice relaxation gener- 
ated by the presence of electrons is much smaller than 
quantum lattice fluctuations. 

In the strong coupling regime (7 = 2t) we have found 
that the optimal states for Nf = 0, 2 are very similar 
to the ground states of oscillators with equilibrium po- 
sitions q » 2Nf-f/ui. This corresponds to a bipolaronic 
ground state in which electrons are trapped by local lat- 
tice distortions and form pairs localized on a single site, 
with one pair for every two sites. Except for a small shift 
in the equilibrium positions, the optimal phonon states 
for a singly occupied site are almost the superposition of 
the optimal phonon states for Nf = 0,2. This can be 
understood as a retardation effect: the Nf = 1 states are 
intermediate states with low probability; electrons do not 
spend enough time in these states for the phonon states 
to adapt. These optimized phonon states for Nf = 1 are 
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very different from the phonon states generated by the 
Lang-Firsov transformation which shifts each osciUator 
equiUbrium position by a quantity 2Nf^/uj. Therefore, 
it is not surprising that the standard strong couphng the- 
ory, which is based on the Lang-Firsov transformation, 
poorly describes the electronic and dynamical properties 
of the Holstein Hamiltonian 

Mean-field theory predicts that the ground state of 
the half-filled Holstein model is a charge-density-wave 
(CDW) state with a dimerized lattice and a gap at the 
Fermi surface for any finite electron-phonon coupling. An 
interesting question is whether the ground state is mod- 
ified when quantum lattice fluctuations are taken into 
account. There is strong evidence that the system is 
dimerized for arbitrary finite coupling at finite phonon 
frequency [§,1). On the other hand, recent results sug- 
gest that the gap is destroyed by quantum lattice fluc- 
tuations in the weak coupling regime I). In the 2-site 
Holstein model it is known that there is a crossover from 
quasi-free electrons to a bipolaron at a finite coupling . 

In Fig. 4 we present several quantities obtained with 
our method as a function of the electron-phonon coupling 
7 for a 6-site system with uj = t. Our results for the 
phonon order parameter nip, defined by Am"^ = {{qe — 

qe+i)'^) ~ 2, where = b^^ + bi, are qualitatively similar 
to the predictions of previous studies Note that 

with this definition, nip = for 7 = 0. We find nip 7^ 
for any finite 7, although nip is smaller than the zero- 
point lattice fluctuations in the weak coupling regime 7 < 
0.8t. For 7 > 1.5t, rUp approaches the strong-couphng 
theory result nip ~ The electronic static staggered 

susceptibility, defined as 



Xf 



(5) 



indicates the existence of a transition around j — t, 
where Xf suddenly increases from the free electron result 
X/ = 1 and approaches the value x/ = 6 representative of 
a perfect CDW order in a 6-site system (Fig. 4). Finally, 
in Fig. 4 we also show the next-nearest neighbor pairing 
correlation (P^P/^2)i where Pi = ci^ci^. This quantity 
has a peak around the value of 7 where the dimeriza- 
tion amplitude nip starts to dominate quantum lattice 
fluctuations and %/ goes up rapidly. 

In the weak coupling regime we have been able to 
study larger systems (with up to 40 sites) with a DMRG 
method using a different approach to handle the phonon 
Hilbert space We have found that despite the peak 
evident in Fig. 4, the pairing correlations {PiP\j^^ de- 
crease as 1/m^, similar to the behavior of free electrons, 
even for a coupling as large as 7 = 0.8i. In the strong cou- 
pling regime, pairing correlations decay exponentially be- 
cause the dimerization opens a gap at the Fermi level [|j . 
Therefore, it is possible that there is a transition from a 
metallic to an insulating CDW phase at a finite electron- 
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FIG. 4. Phonon order parameter mp, staggered static elec- 
tronic susceptibility x/ and next-nearest neighbor pairing cor- 
relations {PiP\j^'^ as a function of the electron-phonon cou- 
pling in a 6-site Holstein lattice with uj = t. 

phonon coupling 

Obviously, a complete understanding of the Holstein 
model properties requires the study of larger systems. 
The use of optimized phonon basis coupled to powerful 
numerical methods, such as DMRG, will enable us to 
perform these calculations. Similarly, the techniques de- 
scribed in this letter could greatly improve our capability 
to perform numerical studies of other problems which in- 
volve an infinite Hilbert space. 
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